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ABSTRACT 

Presented  here  are  the  main  techniques  necessary  to  understand  rotations  in 
three  dimensions,  for  use  with  global  visualisation  and  aerospace  simulations. 
Relevant  techniques  can  be  extremely  difficult  to  find  in  textbooks,  so  some 
useful  examples  are  collected  here  to  highlight  these  techniques. 

The  three  standard  aerospace  coordinate  systems  are  described  and  built 
using  rotations.  The  mathematics  of  rotations  is  described,  using  both  matrices 
and  quaternions.  The  necessary  calculations  are  given  for  analysing  standard 
scenarios  that  involve  the  Global  Positioning  Satellite  system  for  finding  line- 
of-sight  directions  on  Earth,  as  well  as  for  visualising  the  world  from  a  cockpit, 
and  for  converting  to  and  from  the  standard  software  protocol  for  distributed 
interactive  simulation  environments. 

Appendices  then  discuss  combining  rotations,  conversions  with  a  particular 
type  of  Euler  angle  convention,  the  dangers  of  confusing  Euler  angles  with 
incremental  rotations  for  software  writers,  and  finally  there  is  a  short  dicussion 
of  interpolation  of  rotations  in  computing. 
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Using  Rotations  to  Build  Aerospace  Coordinate  Systems 


EXECUTIVE  SUMMARY 

This  paper  presents  the  main  techniques  necessary  to  understand  three-dimensional  ro¬ 
tations.  The  subject  is  not  new,  but  can  be  very  difficult  to  sort  out  and  to  explore  in 
textbooks.  Mathematics  and  physics  texts  that  discuss  the  subject  generally  do  so  only  in 
passing,  and  probably  never  cover  any  standard  software  protocols.  The  main  literature 
on  the  subject  seems  to  be  that  generated  by  the  graphics  computing  community,  since 
programmers  need  to  use  rotations  when  writing  three-dimensional  visualisation  code. 
Nevertheless,  code  ease-of-use  and  ability  to  be  passed  on  means  that  although  there  are 
many  internet  sites  that  talk  about  rotations,  their  information  content  is  often  biased 
towards  whatever  system  or  code  their  authors  have  become  familiar  with. 

We  begin  by  describing  the  three  standard  coordinate  systems  that  are  used  for  sim¬ 
ulation  of  aerospace  scenarios  and  the  Global  Positioning  Satellite  system,  and  how  their 
defining  axes  can  be  built  from  rotations  alone.  We  then  describe  the  mathematics  of 
rotations  using  both  matrices  and  quaternions,  defining  Euler  angles,  and  concentrating 
on  the  important  matrix  (or  equivalently,  quaternion)  that  allows  any  rotation  about  any 
axis  to  be  made. 

As  examples  of  the  techniques,  we  give  the  necessary  calculations  for  dealing  with 
standard  scenarios  involving  the  Global  Positioning  Satellite  system,  for  finding  line-of- 
sight  directions  on  Earth,  for  visualising  the  world  from  a  cockpit,  and  for  converting  to  and 
from  the  standard  software  protocol  for  distributed  interactive  simulation  environments. 

Appendices  then  discuss  combining  rotations,  conversions  with  a  particular  type  of 
Euler  angle  convention,  the  dangers  of  confusing  Euler  angles  with  incremental  rotations 
for  software  writers,  and  finally  there  is  a  short  dicussion  of  interpolation  of  rotations  in 
computing. 
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Glossary,  Conventions,  and  Constants 

a  Latitude  of  a  point  in  the  ECEF  frame. 

l j  Longitude  of  a  point  in  the  ECEF  frame. 

Aircraft’s  own  axes  x,y,z. 

Aircraft’s  own  axes  as  vectors  in  ECEF  x,y,z. 

Earth’s  elliptical  cross-section  This  is  defined  by  the  following  numbers: 

Semi-major:  a  =  6,378,137  nr;  semi-minor:  b  =  6,356,752.3142  m. 

ECEF  Earth-Centred,  Earth-Fixed  Frame.  A  rotating  frame  fixed  to  Earth;  the  global 
frame  of  this  report  within  which  all  computations  are  done. 

ECEF  Cartesian  Axes  X,  Y ,  Z . 

Three  Euler  rotations  act  on  a  set  of  orthonormal  vectors  x0 ,y0,z0  to  give,  in  order, 
then  x2,y2,z2 ,  and  finally  x3,y3,z3. 

GPS  Global  Positioning  Satellite  System. 

Local  Geographic  Frame  A  frame  based  at  our  current  position,  usually  on  or  near 
Earth’s  surface.  Two  common  choices  for  its  axes  are: 

NED  North-East-Down  set  of  axes. 

ENU  East-North-Up  set  of  axes. 

n  A  unit  vector  pointing  along  an  axis  of  rotation.  A  column  when  used  with  matrices, 
and  a  row  when  used  with  quaternions.  This  might  seem  potentially  confusing,  but 
there’s  no  ambiguity  in  practice;  on  the  contrary,  to  continually  write  n  and  rfi  (the 
transpose)  just  to  emphasise  a  column  or  row  would  be  tedious. 

North,  East,  Down  axes  N,E,D. 

North,  East,  Down  axes  as  vectors  in  ECEF  N,E,D. 

Rotation  matrix  Rn{9)  rotates  by  angle  9  around  unit  vector  n  using  the  right-hand 
convention. 

Rotation  quaternion  Qn{9 )  rotates  by  angle  9  around  unit  vector  n  using  the  right- 
hand  convention. 

Position  of  point  A  relative  to  point  B\  i.e. ,  a  vector  pointing  from  B  to  A. 

SLERP  Spherical  Linear  Interpolation.  Used  for  interpolating  orientations  using  quater¬ 
nions. 

WGS-84  World  Geodetic  System  1984  standard  for  maps.  This  standard  defines  the 
oblate  spheroid  typically  used  to  model  Earth’s  shape,  as  used  in  this  report. 
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1  Introduction 


Simulating  air  scenarios  often  involves  the  ability  to  implement  changes  in  orientation. 
The  researcher  might  need  to  visualise  changes  in  orientation,  or  describe  the  departure 
from  flying  “straight  and  level”  for  an  aircraft  flying  at  an  arbitrary  orientation,  undergoing 
various  heading  changes,  pitches,  and  rolls.  Furthermore,  these  calculations  are  dependent 
on  the  aircraft’s  position  on  an  ellipsoidal  Earth,  so  that  the  modeller  has  a  potential 
problem  with  various  rotations  that  need  to  be  performed. 

Add  to  this  the  choice  of  different  ways  to  rotate,  and  it’s  not  surprising  that  rotations 
in  three  dimensions  tend  to  be  seen  as  overly  difficult.  To  complicate  the  situation  further, 
there  are  any  number  of  internet  graphics  software  sites  that  attempt  to  give  information 
and  explanations  about  rotations,  but  these  are  not  always  consistent  in  their  notation, 
and  are  also  not  always  trustworthy.  There  can  be  a  gap  between  the  explanations  of 
those  who  understand  the  mathematics  but  who  might  not  have  any  “on  the  ground” 
programming  experience,  and  those  who  might  not  understand  the  mathematics  but  who 
do  appreciate  the  practical  or  numerical  problems  involved  with  writing  computer  code 
that  performs  rotations.  But  wrongly  interpreting  the  behaviour  of  rotated  objects  can 
lead  to  difficulties,  as  can  be  seen  in  a  later  section. 

This  paper  attempts  to  explain  some  of  the  basics  behind  rotations.  It  is  not  intended  to 
give  lots  of  formulae  for  doing  all  manner  of  rotations,  since  these  are  complicated  enough  to 
be  unusable  if  the  reader  has  no  knowledge  of  just  what  is  happening  underneath.  Instead, 
the  plan  is  to  make  use  of  a  minimum  number  of  techniques  in  order  to  reinforce  the  main 
ideas.  If  some  basic  strategies  applicable  to  solving  common  orientation  scenarios  can  be 
outlined,  then  it’s  hoped  that  the  reader  can  calculate  some  things  from  first  principles. 
If  a  software  package  is  used,  or  even  just  snippets  of  code,  some  knowledge  of  what 
is  really  happening  goes  a  long  way  to  helping  the  user  validate  or  debug  such  code. 
Software  documentation  is  not  above  using  axes  labels  such  as  x,  y,  z  and  X,  Y,  Z,  and 
then  subsequently  mixing  up  the  correct  cases;  so  a  first-principles  knowledge  is  sometimes 
necessary  just  to  read  that  documentation. 


2  Important  Coordinate  Systems 

Understanding  aerospace  simulations  usually  assumes  knowledge  of  three  frames  of  refer¬ 
ence: 


The  Earth-Centred,  Earth-Fixed  Frame  is  global,  attached  to  Earth  itself,  and  al¬ 
ways  uses  both  cartesian  and  polar  coordinates. 

The  Local  Geographic  Frame  is  attached  to  Earth  but  based  at  where  the  aircraft 
currently  is.  It  uses  either  of  two  different  cartesian  coordinate  systems. 


The  Aircraft  Frame  is  fixed  to  the  aircraft  and  uses  three  cartesian  axes,  but  tends  to 
describe  measurements  by  angles  about  these  axes. 
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Figure  1:  The  Earth- Centred,  Earth- Fixed  frame  can  take  a  set  of  cartesian 
axes  X,Y,Z  with  their  origin  at  Earth’s  centre.  The  X-axis  points  to  the 
0°  latitude,  0°  longitude  point.  The  Y -axis  points  to  0°  latitude.  90°  lon¬ 
gitude.  The  Z-axis  points  to  90°  latitude,  along  Earth’s  axis  of  rotation. 
Earth’s  surface  is  taken  to  be  an  oblate  spheroid,  i.e.  having  circular  cross 
section  at  all  latitudes,  and  a  constant  elliptical  cross  section  through  any 
meridian 


There  is  plenty  of  scope  here  for  difficulties  and  obscurity.  In  this  section  we  describe  each 
of  these  frames  and  coordinate  choices  in  turn,  before  discussing  how  to  go  from  one  to 
another. 


2.1  ECEF:  the  Earth-Centred,  Earth-Fixed  Frame 

For  most  situations,  and  certainly  all  in  this  report,  the  Earth-Centred,  Earth-Fixed  Frame 
(ECEF)  is  a  very  useful  stage  on  which  to  play  out  all  scenarios,  from  the  point  of  view 
of  calculating  orientations  and  directions  for  aircraft  with  global  movements.  All  ma¬ 
nipulations  of  numbers  are  done  by  using  cartesian  coordinate  axes  that  are  fixed  in  the 
ECEF. 

The  axes  of  the  ECEF  are  right  handed  and  have  their  origin  at  the  centre  of  Earth, 
being  fixed  to  its  body:  they  rotate  with  it.  They  are  quite  adequate  for  describing 
situations  on  Earth,  but  are  not  so  useful  for  describing  satellite  motion,  such  as  scenarios 
that  involve  the  GPS  satellite  system.  The  reason  is  because  the  ECEF  frame  is  not  quite 
inertial:  because  of  its  rotation,  Newton’s  laws  don’t  quite  hold  in  it  (unless  we  introduce 
complicated  centrifugal  and  Coriolis  forces),  which  means  that  satellite  motion  can  be  very 
complicated  in  the  ECEF’s  coordinates.  For  such  motion,  a  more  encompassing  frame  tied 
to  the  fixed  stars  is  used,  but  we  won’t  need  such  a  one  in  this  report. 

The  ECEF  has  two  common  coordinate  systems:  a  polar-type  “latitude-longitude- 
height”  called  geodetic  coordinates ,  and  the  simpler  three  cartesian  axes  X,Y,Z  that  are 
shown  in  Figure  1.  Global  positions  are  often  given  in  lat-long— height  coordinates,  but 
since  X,  Y,Z  are  far  easier  to  work  with,  we’ll  convert  all  lat-long-height  to  X,  Y,Z.  A  shape 
for  Earth  commonly  used  in  calculations  is  the  oblate  spheroid  specified  by  the  World 
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Geodetic  System  1984  standard  (WGS-84),  which  has  a  circular  cross  section  at  any  given 
latitude,  and  whose  cross  section  through  any  meridian  is  an  ellipse,  having  identical  axes 
lengths  for  all  longitudes.  These  axes  lengths  are,  by  definition, 


Semi-major:  a  =  6,378,137  m 

Semi-minor:  b  =  6,356,752.3142  m.  (2-1) 


These  numbers — especially  b  with  its  extra  decimal  places — might  appear  to  be  specifying 
a  completely  over-optimistic  accuracy;  but  they  are  simply  a  combination  of  measurement 
and  best  fit  to  an  ellipse,  and  so  they  just  define  WGS-84. 

Any  point  has  latitude  a,  longitude  uj,  and  height  h  in  this  report.  (There  is  no 
standard  convention  for  latitude  and  longitude  angles,  so  we  have  used  the  second  letter  of 
each  word  [“a”  and  “o”]  and  written  them  with  Greek  letters.)  The  latitude  and  longitude 
are  shown  in  Figure  2,  while  the  height  of  the  point  is  the  distance  above  the  reference 
spheroid,  i.e.  along  a  line  normal  to  it,  not  along  a  line  extending  from  Earth’s  centre. 

The  reason  the  latitude  is  defined  with  reference  to  the  local  normal  and  not  Earth’s 
centre,  is  because  by  using  the  local  normal,  we  can  be  sure  that  if  two  points  on  the  same 
meridian  have  latitudes  of  say  10°  and  50°,  then  the  normal  line  (i.e.  the  local  vertical)  is 
guaranteed  to  rotate  through  40°  when  passing  from  one  to  the  other.  So  it’s  very  easy  to 
compute  how  the  vertical  changes  over  Earth’s  surface,  and  this  would  not  be  the  case  if 
latitude  was  measured  relative  to  Earth’s  centre. 

Given  lat-long- height  values  for  any  point,  the  corresponding  X,Y,  Z  values  are  given 
in  terms  of  a  and  b  by: 


X 


Y 


Z 


cos  a  cos  uj  , 


cos  a  sin  uj  , 


sin  a . 


(2.2) 


The  inverse  transformation  process  is  lengthier.  Finding  the  longitude  is  the  easiest  part: 


coscu  = 


VX2  +  Y2  ’ 


sin  a;  = 


Y 

VX2  +  Y2  ’ 


(2.3) 


so  that  for  example  in  C  and  Matlab,  uj  can  be  found  by  using  the  single  command 
omega  =  atan2(Y,X).  (Note  that  in  Excel,  the  same  function  atan2  has  its  arguments 
reversed,  so  would  be  atan2(X,Y)  here.) 


Latitude  is  more  difficult,  but  can  be  calculated  iteratively  with  reasonable  ease.  Begin 
with  a  first  estimate 


a  —  tan 


-l 


b2  y/ X 2  +  Y2 


(2.4) 
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where  e.g.  Matlab’s  atan  function  will  suffice  (atan2  is  not  necessary),  and  then  iterate 
to  refine  this  estimate  using 


2  •  2 

a  sin  a 

b2  sin  a  cos  a  +  X 2  +  Y2  sin  a  —  Z  cos  y/ a2  cos2  a  +  b2  sin2  a 

(This  will  fail  if  Z  =  0,  but  then  a  =  0  trivially  anyway.)  This  expression  for  a  converges 
very  quickly;  in  fact  about  five  decimal  places  of  precision  is  obtained  by  only  using  the 
first  estimate  (2.4)  and  not  even  iterating  at  all.  Finally,  after  an  acceptable  estimate  of  a 
has  been  obtained,  the  height  is  calculated: 


Vx2 Ty2 

cos  a 


\/ a 2  cos2  a  +  b2  sin2  a 


(2.6) 


2.2  The  Local  Geographic  Frame 

At  any  point  on  Earth’s  surface,  the  local  geographic  frame  is  defined  by  the  almost  flat 
ground  and  the  vertical  direction,  and  is  relevant  because  it  is  the  basic  reference  the 
aircraft  flies  against,  defining  straight  and  level  flight  and  of  course  the  direction  of  down. 
This  is  an  intuitive  frame  for  any  scenario  whose  total  extent  is  no  more  than  some  tens 
of  kilometres. 

Two  coordinate  systems  are  used:  both  cartesian  with  their  origins  at  the  point  in 
question.  In  the  North-East-Down  or  NED  frame,  shown  in  Figure  2,  the  local  directions 
of  north,  east,  and  down  define  a  right-handed  set  of  three  axes.  An  alternative  choice 
of  cartesian  axes  is  the  local  East-North-Up  set  (ENU).  In  both  axes  sets,  the  north  and 
east  directions  define  a  plane  tangent  to  Earth’s  surface.  At  any  point,  north  points  along 
the  local  meridian,  while  east  points  along  the  local  small  circle  of  constant  latitude. 

The  up/down  or  vertical  direction  is  perpendicular  to  this  plane.  In  general,  the  vertical 
direction  does  not  intersect  Earth’s  centre.  Rather,  it  makes  an  angle  of  the  latitude  with 
the  equatorial  plane,  shown  as  the  angle  a  in  Figure  2. 


2.3  An  Aircraft’s  Own  Frame 

The  aircraft  itself  has  its  own  frame  described  by  a  set  of  axes  at  rest  relative  to  it,  as 
shown  in  Figure  3.  Conventionally  the  x-axis  points  forward  along  the  nose,  the  y- axis 
points  out  along  the  starboard  wing,  and  the  2-axis  points  down.  The  x,  y,  z  axes  of  an 
aircraft  flying  straight  and  level  due  north  will  match  the  local  NED  axes  (i.e.  x  =  north, 
y  =  east,  2  =  down). 


3  The  Basic  Tool:  Vector  Rotations 
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Switching  between  frames  involves  rotating  vectors  about  specified  axes,  and  just  how 
to  rotate  vectors  forms  one  of  the  core  tools  described  in  this  report.  In  this  section  we 
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Figure  2:  A  plane  tangent  to  Earth’s  hypothetical  spheroidal  surface  contains 
the  local  East  and  North  directions,  which  set  those  axes  for  the  North-East- 
Down  frame.  The  Down  axis  is  normal  to  this  plane.  The  latitude  a  of 
the  point  at  the  origin  of  the  NED  frame  is  set  by  the  angle  at  which  the 
Down  axis  intercepts  Earth’s  equatorial  plane.  Because  Earth  is  modelled 
as  an  oblate  spheroid,  the  Down  axis  only  points  toward  its  centre  when  the 
NED  origin  is  on  the  equator  or  at  the  poles.  An  alternative  set  of  axes  is 
East-North-  Up 

present  the  main  formula  for  doing  this.  Surprisingly,  this  matrix  formula  for  a  general 
rotation  in  three  dimensions  can  be  hard  to  find  in  textbooks.  Before  giving  the  general 
expression,  we’ll  start  from  more  basic  ideas. 


3.1  Matrix  Representation  of  an  Orientation 


How  can  we  represent  the  orientation  of  a  body  in  three  dimensions?  We  can  only  do 
it  in  a  relative  way,  by  specifying  how  the  body  has  been  moved  from  some  initial,  or 
base,  position.  Suppose  that  the  body  is  not  translated,  but  its  orientation  is  changed  in 
some  arbitrary  way  from  this  base  position.  Locate  the  body  relative  to  the  unmoving 
basis  vectors,  and  imagine  the  body  taking  a  copy  of  each  basis  vector  along  with  it  as  it 
changes  orientation.  What  are  these  new  vectors?  If  the  old  ones  are 


h- 4 

'o' 

'O' 

0 

III 

1 

hi 

-S£ 

0 

o 

o 

1— 1 

then  suppose  they  are  transformed  into  column  vectors  i' ,  j' ,  k' .  In  that  case,  the  change 
in  orientation  can  be  described  by  the  matrix 

A  =  [»'  j'  k'}  .  (3.2) 


in  the  sense  that  this  multiplies  i  to  give  i' ,  and  similarly  j ,  k.  This  expression  is  useful  in 
that  it  gives  the  matrix  that  transforms  an  initial  orientation  to  a  final  one  immediately, 
without  our  needing  to  work  out  how  the  final  orientation  was  produced. 
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Figure  3:  The  x,y,z  axes  of  the  aircraft’s  own  frame.  These  are  axes  about 
which  the  aircraft  rolls  (x),  pitches  (y),  and  yaws  (z),  where  the  direction  of 
positive  rotation  for  each  is  given  by  the  right  hand  rule 


If  the  body  is  rigid,  then  the  linearity  of  matrix  multiplication  ensures  that  this  matrix 
will  describe  the  orientation  of  the  whole  body.  That  is,  any  other  point  in  the  base- 
positioned  body  that’s  described  by  a  certain  linear  combination  of  the  basis  vectors  i,j ,  k, 
will  be  described  by  the  same  linear  combination  of  the  new  basis  vectors  i',j',k'  in  the 
new  orientation.  So  the  matrix  A  multiplies  any  vector  to  produce  its  orientated  form. 


The  same  result  holds  for  the  more  trivial  two  dimensional  case,  where  we  rotate  all 
vectors  through  9  in  the  xy  plane,  positively  about  the  z-axis.  The  basis  vectors  then 
map  as: 


'1' 

cos  9 

•/ 

'O' 

—  sin# 

0 

¥ 

sin# 

=  *  ,  3  = 

1 

> 

cos  # 

(3.3) 


giving  the  familiar  rotation  matrix 


A  =  [%’  j'] 


cos  9  —  sin  9 
sin  9  cos  9 


(3.4) 


3.2  Euler  Angles 

In  two  dimensions  the  above  transformation  cannot  help  but  be  a  rotation,  but  in  three 
dimensions  it’s  not  explicitly  so;  it  simply  produces  a  final  orientation  given  some  initial 
orientation.  Nevertheless,  it  can  be  shown  to  be  producible  by  just  one  rotation  around 
some  axis,  where  that  axis  needn’t  be  any  of  the  x,  y ,  or  z  axes  (and  in  general  won’t  be). 
This  is  known  as  Euler’s  theorem. 

The  literature  actually  spends  more  time  on  the  fact  that  the  final  orientation  can 
be  produced  by  successive  rotations  around  the  space-fixed  x,  y,  z  axes,  through  what  are 
known  as  Euler  angles.  There  are  many  different  conventions  for  specifying  Euler  angles, 
but  in  any  of  them,  three  angles  are  required  to  describe  a  general  change  in  orientation. 
Some  authors  leave  one  axis  out  and  repeat  another,  so  that  the  orientations  are  performed 
around,  for  example,  the  x-axis,  then  y-axis,  then  x-axis  again  (and  all  similar  choices  of 
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three  distinct  rotations  will  also  work,  such  as  z,  y,  z  etc. — but  not  e.g.  z,  y,  y,  since  this 
is  really  only  two  distinct  rotations).  Some  rotate  around  the  x-axis,  then  the  new  y- axis 
(better  called  y'),  then  the  new  x-axis  (now  called  x'),  although  this  can  be  shown  to 
reduce  to  a  similar  expression  involving  the  same  angles  but  in  a  different  order,  around 
the  space-fixed  x,  y,  x  axes.  The  possibilities  are  many  and  confusing. 

Two  Euler  angles  are  also  used  in  the  field  of  aircraft  tracking,  where  for  example  a 
radar  might  be  used  to  follow  the  motion  of  an  aircraft.  At  any  moment,  the  aircraft’s 
position  can  be  specified  by  azimuthal  and  elevation  angles,  as  if  we  were  using  a  telescope 
with  a  standard  tripod  to  sight  the  aircraft.  When  the  aircraft  is  nearly  overhead,  the 
telescope  becomes  difficult  to  move  easily,  because  the  closer  we  come  to  pointing  it  at 
the  zenith,  the  more  its  azimuthal  motion  (motion  about  the  vertical  axis)  constricts  its 
tube  to  making  smaller  movements.  This  means  that  small  movements  around  the  zenith 
translate  to  large  changes  in  the  azimuthal  angle,  which  can  be  difficult  to  achieve  evenly. 
Mechanically  this  constricted  movement  is  called  gimbal  lock.  In  using  these  two  angles 
to  track  an  aircraft,  rapid  changes  in  azimuth  can  be  difficult  to  model  in  any  numerical 
algorithm,  and  the  term  gimbal  lock  has  also  come  to  be  applied  to  numerical  difficulties. 
The  result  is  that  it  has  given  Euler  angles  a  bad  name.  But  the  problem  is  purely  one  of 
numerical  stability;  there’s  nothing  wrong  with  using  the  two  Euler  angles  in  principle. 

This  constriction  of  motion  about  one  axis  also  causes  mechanical  problems  with  some 
gyroscopes.  To  picture  why  this  might  be,  imagine  a  very  simple  (and  not  very  useful) 
example  of  an  inertial  navigation  system:  a  horizontal  spinning  flywheel  carried  on  a  yoke 
attached  to  hinges  at  the  end  of  each  wingtip,  and  held  underneath  the  aircraft  while  it 
flies  straight  and  level.  The  flywheel’s  axis  determines  the  up/down  direction  used  by 
the  autopilot  to  help  keep  the  aircraft  on  course.  The  aircraft  can  pitch  up  and  down 
and  yaw  in  the  local  north-east  plane  without  affecting  the  flywheel’s  motion.  But  if  the 
human  pilot  tries  to  roll,  the  flywheel  will  be  forced  out  of  its  spin  plane.  In  such  an  event, 
the  inertial  navigation  system  will  be  fooled  into  believing  that  the  aircraft  has  begun 
to  develop  some  serious  instabilities  which  it  will  try  to  correct,  with  perhaps  disastrous 
consequences.  While  this  is  a  big  simplification  of  a  real  system,  it  does  encapsulate  the 
problem  that  can  happen  with  a  real  gyroscope  when  the  plane’s  attitude  reaches  some 
extreme  position. 

One  way  of  avoiding  numerical  instabilities  in  tracking  algorithms  is  to  use  Euler’s 
theorem,  which  does  away  with  having  to  use  three  angles  around  the  three  fixed  axes.  So 
it  is  that  our  matrix  A  in  (3.2)  describes  this  rotation.  In  principle  it  should  be  sufficient 
for  all  purposes,  although  in  practice  its  nine  entries  can  require  a  powerful  ability  for 
number  crunching,  especially  when  used  to  describe  motions  where  the  orientation  angles 
are  changing  with  time. 


3.3  The  Workhorse:  Rotating  in  Three  Dimensions 

The  matrix  A  in  (3.2)  describes  the  rotation  that  changes  a  body’s  orientation,  provided 
we  know  how  the  basis  vectors  change.  But  usually  we  are  not  given  the  final  position 
of  those  vectors.  In  general,  we  need  to  be  able  to  rotate  an  arbitrary  vector  around  an 
arbitrary  axis,  and  this  is  a  slightly  different  task  to  be  done. 
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Figure  4:  The  matrix  Rn(9 )  rotates  x  in  a  right-handed  sense  around  the 
unit  vector  n  by  an  angle  6  to  produce  x' .  That  is,  x'  =  Rn(9 )  x 


This  rotation  can  be  done  by  multiplying  the  vector  by  a  matrix  Rn{9)  that  we  write 
down  in  this  section.  Rn{9)  specifies  a  right-handed  rotation  through  some  angle  9  about 
an  axis  aligned  with  a  unit  vector  n.  (Requiring  n  to  have  unit  length  makes  the  expres¬ 
sions  simpler.)  It  rotates  the  vector  x  to  give  x' ,  as  shown  in  Figure  4: 

x'  =  Rn{9)  x  .  (3.5) 

Write  the  unit  vector  that  points  along  the  axis  of  rotation  as 


n  = 


(3.6) 


where  n  can  point  in  either  direction  along  the  axis,  but  changing  the  choice  of  that  direc¬ 
tion  will  change  the  sign  required  for  the  rotation  angle  9,  since  the  rotation  matrix  Rn(9) 
obeys  the  right  hand  rule  for  rotation.  The  rotation  matrix  can  be  written  in  the  following 
way,  which  is  the  main  equation  of  this  report: 


n\  nln2  n1n3 

0 

-n3 

n2~ 

Rn{9)  =  (1  —  COS  9) 

n2n1  n2  u2n3 

+  cos  9  I3  +  sin  9 

n3 

0 

~ni 

n3nl  n3n2  n3 

-n2 

n\ 

0_ 

=  (1  —  cos  9)  nn1  +  cos  9  I3  +  sin  9  nx  ,  (3.7) 

where  n t  is  the  transpose  of  n,  I3  is  the  3x3  identity  matrix,  and 


n 


X 


0  — n3  n2 

n3  0  -rq  , 
— n2  rq  0 


(3.8) 


so-named  because 


nxx  =  n  x  x  , 


(3.9) 


(with  the  n  written  nonbold  to  emphasise  its  matrix  form).  The  rotation  matrix  Rn(9) 
is  called  orthogonal ,  by  which  is  meant  RR!  =  RfR  =  1,  and  as  well  its  determinant  is 
always  1.  Equations  (3.5)  and  (3.7)  will  be  used  repeatedly  in  this  report  to  break  more 
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complicated  procedures  up  into  single  rotations  that  are  easy  to  visualise  and  easy  to  code. 
They  are,  in  fact,  the  central  equations  of  the  whole  of  rotation  theory. 


Example  of  the  use  of  (3.7):  Rotate  the  vector  (2,0,0)  by  90°  about  the  y-axis. 
What  vector  results?  The  required  rotation  matrix  is  Ry( 90°)  (where  by  the  subscript  y 
is  meant  n  =  (0, 1,  0)*).  Equation  (3.7)  gives  it  as 


Ry(  90°)  =  nof  +  nx 


'O' 

1 

o 

0 

1' 

1 

o 

0 

1" 

1 

[0  1  0]  + 

0 

0 

0 

= 

0 

1 

0 

o 

-1 

0 

CD 

-1 

0 

0 

(3.10) 


(Actually,  in  this  simple  example  we  can  calculate  72^(90°)  alternatively  using  (3.2).  Simply 
note  that  the  basis  vectors  rotate  as 


1 

0 

0 

0 

0 

1 

0 

— > 

0 

5 

1 

— > 

1 

5 

0 

— > 

0 

0 

-1 

0 

0 

1 

0 

so  that  (3.2)  yields  Ry( 90°)  trivially.)  The  required  rotated  vector  is  then 


72,(90°) 


'cm  o 

1 _ 

_ 

t  o 

0 

1' 

'2' 

o' 

1 

0 

0 

= 

0 

0 

0 

i  o 

-2 

(3.11) 


(3.12) 


as  expected. 


An  example  of  combining  two  rotations  is  given  in  Appendix  A.  As  stated  by  Euler’s 
theorem,  the  effect  of  these  is  just  one  rotation,  although  the  values  of  the  resulting 
equivalent  axis  and  angle  turned  through  are  not  obvious  at  all  from  the  original  two  axes 
and  angles. 


The  rows  of  a  rotation  matrix  will  gradually  lose  their  mutual  orthonormality  because 
of  rounding  errors  under  repeated  iterations  in  a  computer  programme.  This  wandering 
can  be  kept  in  check  by  periodically  re-orthogonalising,  accomplished  by 

72  — >  72  (  72*72) _1/2.  (3.13) 


Rotation  Order  and  its  Possible  Blind  Application 

Since  matrix  multiplication  is  not  commutative,  two  rotations  are  also  not  in  general  com¬ 
mutative:  the  result  depends  on  the  order  in  which  they  are  done.  However,  infinitesimal 
rotations  are  commutative.  This  fact  can  be  used  to  our  advantage  when  writing  com¬ 
puter  code  that  takes  inputs  from  a  joystick,  in  order  to  rotate  a  scene  (say,  in  a  flight 
simulator).  The  pilot  might  induce  both  a  pitch  and  a  roll,  but  only  by  applying  each  of 
these  in  tiny  increments  will  the  software  successfully  reproduce  the  effect  that  the  pilot 
wants. 

Nevertheless,  it’s  possible  to  go  to  an  extreme  of  demanding  a  set  rotation  order  that 
does  not  mimic  what  the  user  of  a  software  package  wants,  but  that  does  do  something 
highly  misleading.  In  fact,  a  search  through  the  internet  shows  that  this  mistake  seems  to 
be  a  classic  of  computer  graphics  programming,  and  has  caused  some  in  that  community  to 
distrust  the  mathematics  of  rotation.  An  example  of  this  incorrect  application  of  rotations 
is  discussed  in  Appendix  C. 
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3.4  Rotating  using  Quaternions 

Although  the  rotation  matrix  Rn{9)  has  nine  entries,  it’s  really  only  built  from  three  num¬ 
bers:  the  angle  9  turned  through,  and  any  two  components  of  n;  the  third  component  of  n 
is  then  implied,  since  n  has  unit  length.  Using  Rn(9)  can  be  computationally  inefficient, 
since  not  only  do  all  nine  components  need  to  be  manipulated,  but  if  the  matrix  is  ap¬ 
plied  within  a  loop,  perhaps  thousands  of  times,  numerical  inaccuracies  can  degrade  its 
orthogonality.  That  is,  its  rows  (or  equivalently  columns)  will  slowly  lose  their  mutual 
orthonormality.  It  can  certainly  be  periodically  re-orthogonalised  after  an  appropriate 
number  of  iterations,  as  shown  in  (3.13).  Nevertheless,  the  question  arises  as  to  whether 
any  way  exists  of  rendering  the  matrix  down  to  its  basic  three  numbers,  and  perhaps 
eliminating  some  of  the  numerical  complexity  in  the  process. 

We  are  always  free  to  construct  any  quantity  from  those  three  numbers.  For  instance  we 
could  specify  a  rotation  through  9  around  (n1,  n2,  n3)  by  a  new  object  written  as  ( 9 ,  n1,n2). 
This  is  very  concise,  but  the  catch  is  that  we  would  also  have  to  define  how  this  object 
acts  on  a  vector  to  rotate  it.  The  result  is  not  anything  simple,  so  we  gain  nothing  by 
doing  this. 

But  all  is  not  lost.  If  we  allow  a  little  more  complexity,  and  construct  a  new  ob¬ 
ject  Qn(9 )  from  9 ,  n1,n2,n3  along  with  a  sine  and  cosine,  then  this  turns  out  to  be  simple 
enough  to  offset  the  more  complicated  way  it  has  to  interact  with  vectors.  Here  n  is  taken 
to  be  the  unit  row  vector  (n1,n2,n3): 

Qn{9)  =  (cos^,nsinQ  .  (3.14) 

This  new  object  Qn{9)  is  called  a  rotation  quaternion.  More  general  quaternions  (not  for 
rotating)  are  composed  of  a  number  and  a  vector:  (o0,a).  Two  quaternions  multiply  in 
the  following  way: 


(a0,  a)  (b0,b)  =  (a0b0  —  a-b,  a0b  +  b0a  +  axb) ,  (3.15) 

so  that  the  squared  length  of  a  quaternion  is  defined  to  be  (in  analogy  to  the  length  of  a 
vector): 

|(°o,  a)|2  =  ao  +  la|2  •  (3.16) 

We  can  see  straight  away  that  a  rotation  quaternion  Qn{9)  has  unit  length,  which  is  the 
equivalent  property  of  rotation  quaternions  as  orthogonality  and  a  unit  determinant  is  of 
rotation  matrices. 

Using  rotation  quaternions,  a  row  vector  x  can  be  rotated  through  an  angle  9  around 
a  unit  vector  n  by  applying  a  double  multiplication: 

(O,a')  =  Qn(0)(O,aOQ-„(0),  (3.17) 

where  quaternion  multiplication  is  associative,  so  either  product  can  be  done  first.  As 
an  example,  in  (3.12)  we  used  a  matrix  to  rotate  (2,0,0)  by  90°  about  the  y-axis  to 
give  (0,0,  —2).  Contrast  the  matrix  approach  with  the  quaternion  approach: 

("’O  =  (^>'1'0))  <0-2-°-0)  (^>-1-0)) 
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=  (1,0, 1,0)  (0,1, 0,0)  (1,0, -1,0) 

=  (0,0, 0,-2),  (3.18) 

so  that  the  result  is  (0, 0,  —2)  as  expected. 

Combining  two  rotations  using  quaternions  is  straightforward.  Just  as  a  matrix  ro¬ 
tation  R\  followed  by  7?  2  is  equivalent  to  a  single  matrix  rotation  R2  R±,  a  quaternion 
rotation  Q\  followed  by  Q2  is  equivalent  to  a  single  quaternion  rotation  Q 2  Q\. 

Finally,  just  as  a  rotation  matrix  will  slowly  lose  orthogonality  under  repeated  iterations 
in  a  software  routine,  so  a  rotation  quaternion’s  length  will  slowly  wander  from  one  in  the 
same  circumstances.  It  can  be  periodically  reset  by  dividing  each  of  the  quaternion’s  four 
elements  by  the  length  calculated  from  (3.16),  which  is  a  very  much  simpler  task  than 
re-orthogonalising  a  matrix  in  (3.13). 


Matrix  Cost  versus  Quaternion  Cost 

How  expensive  is  using  matrices  as  opposed  to  using  quaternions?  While  (3.18)  might 
look  succinct,  it  does  hide  some  effort,  especially  in  the  two  cross  products  that  have  been 
done.  In  Table  1  we  give  simple  counts  of  operations  required  to  rotate  a  vector,  if  done 
with  pen  and  paper.  No  allowance  for  numerical  optimisation  is  made  since  that  forms  a 
field  in  itself.  It  can  be  seen  that  quaternions  are  easier  to  build,  while  matrices  are  easier 
to  use;  thus,  depending  on  the  application,  it  can  be  useful  to  convert  between  the  two 
subject  to  how  much  of  each  process  needs  to  be  done. 


Table  1:  First-principles  computation  costs  of  matrices  and  quaternions,  for 
the  build  and  one  use  of  each.  The  numbers  shown  are  not  optimised 


Matrix 

Quaternion 

Build 

using  (3.7)  and  (3.14): 

2  trigonometric  functions, 
10  additions, 

21  multiplications. 

2  trigonometric  functions, 
no  additions, 

4  multiplications. 

One  rotation 

using  (3.5)  and  (3.15,  3.17): 

6  additions, 

9  multiplications. 

24  additions, 

30  multiplications. 

Multiplication 

27  multiplications, 

18  additions. 

16  multiplications, 

12  additions. 

The  benefit  of  using  quaternions  for  iterative  calculations  is  that  the  resetting  of  a 
quaternion’s  length  to  1  requires  a  simple  division  by  its  length,  calculated  from  (3.16). 
On  the  other  hand,  re-orthogonalising  a  matrix  from  (3.13)  is  a  much  more  complicated 
affair  that  involves  the  lengthy  procedure  of  matrix  diagonalisation.  Although  there  is 
room  for  optimisation  in  this  process,  it  cannot  compete  with  quaternion  normalisation  for 
simplicity.  How  often  such  a  resetting  or  re-orthogonalising  needs  to  be  done  is  presumably 
a  function  of  the  numerical  complexity  of  the  problem  being  tackled. 
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In  fact,  it  might  be  thought  that  rather  than  use  the  expensive  (3.13)  to  re-orthogonalise 
the  matrix,  that  a  better  approach  would  be  to  convert  the  matrix  to  a  quaternion,  nor¬ 
malise  that,  and  then  convert  back  to  a  matrix.  This  can  be  done  using  the  analysis  of 
Appendix  A,  by  applying  (Al)  to  produce  the  angle  and  axis  from  the  matrix,  that  are 
then  used  to  build  the  quaternion  via  (3.14),  which  is  then  normalised  and  converted  back 
to  a  matrix  using  (3.7).  But  whether  this  is  really  viable  might  depend  on  the  numerical 
accuracy  of  the  software.  As  it  is,  (Al)  only  samples  the  matrix’s  three  diagonal  elements 
to  calculate  the  angle  9,  perhaps  producing  too  inaccurate  a  result  to  be  used  in  building 
the  quaternion.  The  quaternion  could  be  produced  by  using  all  of  the  matrix  elements;  but 
it’s  not  clear  what  work  has  been  done  in  this  area  for  doing  this  efficiently  and  accurately. 


4  Rotating  Axes  to  Change  Coordinates 

The  building  block  that  allows  all  manner  of  aerospace  re-orientations  to  be  done,  for 
various  pitches,  rolls,  yaws,  changes  of  latitude/longitude,  and  local  geographic  frames,  is 
the  act  of  rotating  each  of  the  vectors  that  represent  the  axes  of  some  frame. 

There  are  two  basic  tasks  to  consider  in  this  type  of  aerospace  calculation.  First,  given 
a  place  on  Earth  (usually  in  lat-long-height  coordinates),  we  need  to  construct  the  local 
geographic  frame’s  axes,  such  as  NED.  Second,  an  aircraft  can  be  introduced  into  this 
NED  frame,  and  its  orientation  found  relative  to  that  frame. 


4.1  Constructing  NED  Axes  for  the  Local  Geographic  Frame 

Given  the  lat-long-height  of  a  place  on  or  near  Earth,  the  first  thing  we  wish  to  do  is  find 
the  local  NED  axes.  The  technique  used  is  the  primary  one  of  this  report.  We  build  an 
initial  set  of  NED  axes  in  a  place  where  it’s  simple  to  determine  them,  and  then  we  rotate 
these  around  to  the  required  place. 

The  simplest  place  to  construct  an  initial  set  of  NED  axes  is  on  the  junction  of  the 
Equator  and  the  prime  meridian,  since  this  has  a  latitude  and  longitude  of  a  =  uj  =  0°. 
Represent  each  axis  by  a  unit  vector  in  the  ECEF  frame,  which  all  calculations  are  being 
done  within: 


"O' 

'O' 

'-1' 

No  = 

0 

,  Eo  = 

1 

j  Do  = 

0 

1 

0 

0 

The  scenario  is  shown  in  Figure  5.  Now  rotate  each  of  the  No,  E q,  D q  vectors  in  two  steps. 
The  first  rotation  is  about  the  IV  o  vector  by  the  longitude  u>,  taking  care  to  realise  that 
the  rotation  matrix  will  obey  the  right  hand  rule.  This  rotation  creates  an  intermediate 
triplet  of  vectors,  taking  Nq-*Ni,  Eq—*E\,  Dq—>D\.  (Of  course,  IV i  =  N o,  but  for 
clarity  we  keep  each  triplet  together  notationally.) 

Now  rotate  each  of  the  intermediate  set  by  the  latitude  a  about  —E i;  again  we’ve  been 
careful  to  note  the  use  of  the  right  hand  rule,  so  need  to  specify  the  axis  of  rotation  as  —E\. 
(We  can  certainly  use  +E\,  but  will  have  to  rotate  by  —a  about  this.  The  rotation  matrix 
or  quaternion  will  be  unchanged.)  Because  the  original  No,  Eo,  Do  vectors  each  have 
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Figure  5:  Constructing  local  north-east-down  axes  N,E,D  at  a  given  point 
of  known  latitude  and  longitude.  Start  with  the  three  vectors  Nq,  Eq,  Dq. 
Construct  the  intermediate  set  N±,Ei,Di  by  rotating  the  original  set 
through  the  longitude  u)  about  Nq.  Then  rotate  N i,  E\,D\  about  —E\  by  the 
latitude  a,  according  to  the  right  hand  rule,  to  give  the  final  set  N,  E,  D.  Or 
equivalently,  rotate  the  intermediate  set  about  E\  through  minus  the  latitude. 
The  resulting  vectors  are  the  local  north-east-down  axes 


unit  length,  they  automatically  satisfy  the  requirement  for  the  axis  vector  n  to  be  of  unit 
length  in  (3.7).  Notice  that  the  clever  way  that  latitude  has  been  defined  means  we  needn’t 
worry  that  Earth’s  shape  is  not  spherical;  two  places  separated  by  say  50°  latitude  on  the 
same  meridian  will  certainly  have  a  50°  angle  between  their  North  vectors,  and  a  50°  angle 
between  their  Down  vectors. 

The  sequence  of  steps  is: 


N\  =  RNq  (u)  N0  =  N o 

El  =  RNo  (uj)  Eq 

Di  =  RNo(uj)Dq 

N  =  R  Ei  (a)  N i 

E  =  R  E  (a)  Ei  =  Ei 

D  =  R  Ei  (a)  Di  .  (4.2) 

The  resulting  vectors  N,  E,  D  are  the  local  north-east-down  axes,  and  can  be  used  for 
further  calculations.  In  the  interest  of  pedagogy,  we  have  made  the  steps  as  alike  as 
possible,  but  in  practice  they  can  be  heavily  reduced  to 


E  =  RNo(u)Eo 
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AT  =  R_E(a)N  o 

D  =  NxE.  (4.3) 

This  reduction  is  purely  a  question  of  pedagogy  versus  computational  speed.  Equa¬ 
tions  (4.2)  allow  for  a  step  by  step  follow-through  of  the  process,  which  is  useful  for 
writing  transparent  computer  code  or  debugging  in  more  complicated  situations  where 
more  than  two  rotations  need  to  be  performed.  In  this  two-rotation  case,  it’s  certainly  not 
difficult  to  use  (4.3)  without  having  to  think  along  the  lines  of  (4.2)  at  all.  But  situations 
requiring  more  than  two  rotations  will  not  be  so  easily  abbreviated. 

The  following  example  draws  many  of  these  ideas  together. 

Example:  If  we  are  in  Adelaide  and  Earth  is  transparent,  what  is  the  compass  bearing  of 
Brussels  if  we  are  looking  straight  at  it  through  Earth? 

We’ll  calculate  ECEF  position  vectors  of  both  Adelaide  and  Brussels,  subtracting  one 
from  the  other  to  find  the  position  vector  of  Brussels  relative  to  Adelaide,  and  then  express 
this  in  the  local  NED  axes  at  Adelaide  using  dot  products.  It’s  helpful  to  write  the  vectors 
in  the  following  way,  where  A  =  Adelaide,  B  =  Brussels,  and  C  =  some  other  useful  point. 
The  vector  giving  the  position  of  B  relative  to  A  can  be  written  rB<_A,  and  the  following 
are  all  true  in  general,  as  well  as  being  useful  for  writing  down  vector  relationships  without 
needing  to  draw  pictures: 


rB^A 

II 

1 

T 

to 

rB<—A 

—  rB^C  rA^C 

rB^A 

=  rB^C  +  rC<-A 

(4.4) 

The  cities’  positions  are: 

Adelaide:  latitude  a  =  —34.9°,  longitude  u  =  138.5°, 

Brussels:  latitude  a  =  50.8°,  longitude  u  =  4.3°. 

Use  (2.2)  to  find  the  position  of  each  city  relative  to  Earth’s  centre  (which  is  exactly  what 
the  ECEF  coordinates  specify):  set  C  =  Earth’s  centre. 


'-3.92' 

'4.03' 

3.47 

x  106  metres,  rB^c  = 

0.30 

-3.63 

4.92 

x  106  metres. 


The  position  of  Brussels  relative  to  Adelaide  is  then 


'  B< 


-A  ~  rB<—C  rA<—C 


7.95 

-3.17 

8.55 


x  106  metres. 


(4.5) 


(4.6) 


This  vector  is  still  an  ECEF  vector,  although  we  can  consider  it  to  start  at  Adelaide  and 
point  to  Brussels.  To  calculate  the  bearing  of  Brussels  as  seen  from  Adelaide,  we  need 
to  get  the  local  north  and  east  components  of  To  do  this,  we’ll  first  calculate 
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Adelaide’s  local  north,  east,  and  down  axis  vectors,  labelled  N,E,D  in  Figure  5,  all  in 
ECEF  coordinates,  using  (4.2).  For  this  we  need  RNq(uj),  constructed  from  (3.7): 


'O' 

'-0.74896 

-0.66262 

O' 

AT0  = 

0 

,  w  =  138.5°, 

R  (u))  ~ 

JV0  v  > 

0.66262 

-0.74896 

0 

1 

0 

0 

1 

This  RN(j  (w)  then  acts  in  (4.2)  to  produce  N\,  E\,  D\.  We  then  use  —E\  to  construct  the 
next  rotation  matrix  R_  (a),  which  finally  produces  the  N,  E,  D.  The  steps  are  easy  to 
programme  and  are  omitted;  the  resulting  axes  are 


'-0.429' 

'-0.663' 

0.614" 

0.379 

,  E  = 

-0.749 

,  D  = 

-0.543 

0.820 

0 

0.572 

The  bearing  of  Brussels  is  calculated  by  projecting  rB_  4  into  the  local  north-east 
plane.  The  projection  can  be  visualised  as  an  arrow  drawn  on  the  ground,  so  that  its 
angle  from  north  is  the  required  bearing.  The  projection  in  the  local  north-east  plane  is 
simply  given  by  the  components  of  rB^A  on  the  local  north  and  east  axes.  These  are  dot 
products: 

North  component  =  rB^A-N  =  2.4035  x  106  metres, 

East  component  =  rB^  4-  E  =  —2.8958  x  106  metres.  (4.9) 

(Note  that  these  equations  assume  the  axis  vectors  N,  E  have  unit  length;  another  reason 
to  prefer  unit  length  axis  vectors  from  the  start.)  The  units  here  are  immaterial;  the  only 
thing  that  matters  is  the  ratio  of  the  components,  so  that  Brussels’  bearing  is 

360°  -  tan”1  ~  310°  ,  (4.10) 

or  roughly  north-west.  Apart  from  the  initial  change  of  variables  in  (2.2),  this  whole 
calculation  has  been  accomplished  using  just  rotations  and  dot  products.  We  could  of 
course  have  streamlined  it  using  (4.3)  as  well  as  omitting  the  calculation  of  the  local  Down 
vector  D,  but  the  point  here  is  pedagogy,  not  efficiency. 


4.2  The  World  as  Seen  by  a  Pilot 

Sometimes  we  wish  to  know  the  appearance  of  a  scenario  as  seen  by  the  pilot.  For  example 
if  an  aircraft  is  flying  straight  and  level,  and  then  executes  a  series  of  yaws,  pitches,  and 
rolls,  where  will  the  compass  directions  lie,  and  where  will  “up/down”  be?  If  an  aircraft  is 
alongside  but  200  m  higher,  where  will  this  aircraft  be  seen  to  lie  from  the  confines  of  the 
cockpit? 

These  questions  are  answered  by  rotating  vectors.  To  determine  where  the  compass 
directions  are,  we  need  to  know  where  the  aircraft  is  headed  before  it  turns.  These  direc¬ 
tions  are  only  really  well  defined  near  Earth’s  surface,  and  in  this  case,  each  of  them 
can  be  considered  as  vectors  emanating  from  the  cockpit  itself,  described  in  the  air¬ 
craft’s  frame.  For  example,  Fig.  3  shows  that  if  the  aircraft  is  flying  straight  and  level 
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north-east,  then  the  north  direction  vector  will  have  coordinates  in  the  aircraft’s  frame 
of  (. x,y,z )  =  (l/\/2,  — 1/a/2, 0). 

Example:  This  involves  large  distances  just  to  ensure  the  method  is  well  understood.  Sup¬ 
pose  we  are  flying  north-east  over  Adelaide,  pitched  up  at  20°  at  an  altitude  of  30,000  m. 
We  sight  an  aircraft  flying  over  Sydney  at  the  same  altitude.  In  which  direction  in  the 
cockpit  must  we  look  to  see  that  aircraft? 

As  usual,  work  in  the  ECEF.  Call  the  other  aircraft  “them”,  and  construct  a  vector 
pointing  from  our  location  to  them.  This  is  the  vector  of  them  relative  to  us,  or  rthem^us- 
We  require  the  components  of  this  vector  in  the  aircraft  frame.  These  components  are 
given  by  dotting  rthem^ua  with  the  three  vectors  describing  each  of  the  aircraft’s  axes. 
The  axes  are  x,  y,  z  in  Fig.  3,  and  we’ll  call  the  vectors  along  each  one  x,  y ,  2,  respectively. 
So  we  need  to  find 

^thenu — us ' X  ’  ^hem^us^’  ^em^us’*'  (4-U) 

The  vector  rthem<_us  is  found  simply  by  referring  everything  to  the  centre  C  of  Earth, 
which  implies  using  the  ECEF: 

^them< — us  =  rthem^C  “  Tus^C  '  (4'12) 

The  vector  rth  <_c  *s  their  position  in  the  ECEF,  and  TiS<— c  is  our  °wn  position  in 
the  ECEF.  These  positions  are  given  by  applying  (2.2)  to  Adelaide’s  and  Sydney’s  latitude 
and  longitude,  together  with  the  given  heights  both  of  30,000  m.  The  two  cities  lie  at: 

Adelaide:  latitude  a  =  —34.9°,  longitude  uj  =  138.5°, 

Sydney:  latitude  a  =  —33.9°,  longitude  u  =  151.2°. 


Equation  (2.2)  together  with  the  appropriate  heights  then  gives  approximately: 

x  106  metres.  (4-13) 


^therm — C 

'-4.67 

2.57 

x  106  metres, 

^us<— C  = 

'-3.94' 

3.49 

-3.55 

-3.65 

Hence 


rp  —  rp  -  p 

them'* — us  them-* — C  us-* — C 


-7.25 

-9.21 

0.92 


x  10°  metres. 


(4.14) 


Now  we  need  to  find  x ,  y ,  z  at  our  location.  First  we  find  the  local  north,  east,  and 
down  vectors  IV,  E,  D  by  (4.2)  or  (4.3),  where  the  heights  play  no  part  in  the  calculation. 
These  were  given  in  (4.8).  These  vectors  must  be  rotated  to  produce  our  own  aircraft 
axes.  First,  place  our  aircraft  in  the  local  NED  frame  flying  north,  straight  and  level.  Its 
aircraft  axes  then  match  the  local  NED  axes: 


x0  =  N,  y0  =  E ,  z0  =  D.  (4.15) 

Next  rotate  the  aircraft  axes  around  z0  by  45°,  and  then  rotate  20°  around  the  new  y- 
axis  r/j  ,  taking  care  to  get  the  angle  signs  right  by  way  of  the  right  hand  rule  for  rotation. 
The  sequence  is: 


xi  —  RZq  (45  )  x0 
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Vi  —  -RZo(45°)y0 

zi  =  -^z0(45°)z0  =  z0 

x  =  RVi  (20°)  x1 
y  =  RVl  (20°)  yx  =  2/1 
z  =  Ry1  (20  )  Z\  ■ 

The  final  vectors  are  approximately 


'-0.94' 

'-0.17' 

0.31' 

X  = 

-0.06 

,  y  = 

-0.80 

,  *  = 

-0.60 

0.35 

-0.58 

0.74 

(4.16) 


(4.17) 


Dotting  these  with  rt}iem,_lls  hr  (4.11)  gives  the  set  of  components  we  require,  so  that  the 
other  aircraft  turns  out  to  have  coordinates  (765, 802,  393)  km  in  the  frame  of  our  aircraft. 
(Note  that  these  numbers  result  from  keeping  more  significant  figures  than  are  shown  in 
the  numbers  above.)  As  a  simple  check,  the  length  of  this  vector  should  be  the  distance 
of  the  aircraft: 


distance  of  aircraft  =  \/7652  +  8022  +  3932  km  ~  1176  km,  (4-18) 

which  is  reasonable.  To  actually  sight  the  aircraft,  we  first  need  to  turn  our  head  to  our 
right  by  approximately  tan-1  |||,  or  46°.  Again  this  makes  sense,  since  the  aircraft  is 
approximately  due  east  and  we  are  already  flying  north-east.  After  having  turned  our 
head  through  46°,  we  need  to  turn  it  downwards  by  tan-1  ^==p,  or  about  20°. 

4.3  ECEF  Coordinates  and  Heading-Pitch-Roll  Conversion 

Typically,  simulation  software  will  store  object  locations  in  lat-long-height,  as  well  as 
orientations  in  a  heading  pitch-roll  format:  three  Euler  angles  with  respect  to  its  local 
NED  frame.  “Heading”  here  means  where  the  aircraft’s  x-axis  points,  not  the  direction  in 
which  it’s  actually  travelling  over  the  ground,  which  is  called  its  track.  Note  that  heading 
is  not  the  same  as  yaw.  Yaw  is  the  angle  between  an  aircraft’s  heading  (into  wind,  not  over 
ground),  and  the  direction  in  which  its  nose  points.  An  aircraft  usually  flies  with  no  yaw, 
meaning  that  it’s  pointing  exactly  into  the  oncoming  air  flow;  the  air  flow  is  approaching 
directly  over  its  nose.  We,  together  with  software  packages,  are  really  assuming  the  aircraft 
is  not  yawed,  which  is  a  very  normal  way  to  fly.  Only  some  extreme  fliers  will  hold  a  yaw 
after  finishing  their  turn. 

Despite  lat-long-height  and  heading-pitch-roll  formats,  an  IEEE  standard  specifies 
that  in  order  to  communicate  location-orientation  between  different  machines  in  a  dis¬ 
tributed  interactive  simulation  (DIS)  environment,  XYZ  coordinates  be  used,  along  with 
a  different  set  of  Euler  rotations.  The  DIS  standard  says  that  three  Euler  rotations  are 
implemented  by  first  rotating  around  the  2-axis,  then  around  the  new  y-axis,  and  finally 
around  the  newest  x-axis.  This  order  is  useful  because  it  corresponds  to  physical  situations; 
for  example,  ones  in  which  an  aircraft  might  manoeuvre  by  changing  its  heading  (about  its 
2-axis),  then  pitching  (about  its  new  y-axis),  and  finally  rolling  (about  its  newest  x-axis). 


17 


DSTO-TN-0640 


Proprietary  software  might  well  include  routines  to  convert  to  and  from  these  coordinates, 
but  this  software  tends  to  be  provided  without  source  code  and  so  generally  is  not  portable. 
To  aid  in  developing  routines  that  use  the  DIS  standard,  this  section  describes  how  such 
conversions  can  be  done. 

The  first  set  of  DIS  coordinates  is  the  aircraft’s  location  in  the  XYZ  coordinates  of  the 
ECEF  frame.  The  second  set  of  coordinates  is  the  set  of  three  Euler  angles  that  specify 
the  aircraft’s  orientation  by  starting  with  its  own  axes  coincident  with  the  ECEF’s  XYZ 
axes,  and  then  rotating.  The  three  rotations  are  by  angle  ip  around  Z,  then  by  9  around 
the  rotated  Y,  and  finally  by  <p  around  the  latest  rotated  X . 

After  distribution,  these  numbers  (A,  Y,  Z ,  ip ,  0 ,  (f>)  might  need  to  be  converted  back  to 
a  lat  long-height  and  an  orientation  with  respect  to  the  local  NED  axes. 

Converting  the  aircraft’s  location  to  and  fro  between  XYZ  and  lat-long-height  was 
done  in  (2.2)— (2.6).  The  next  task  we’ll  address  is  calculating  the  Euler  angles  ip,  9,  (p.  As 
we’ll  see  in  (4.30)  ahead,  it’s  necessary  to  solve  the  problem  for  a  general  set  of  starting 
vectors,  as  opposed  to  X,Y,Z ,  so  we’ll  restate  it  using  more  general  vectors.  Suppose 
we  have  an  orthonormal  set  of  vectors  (i.e.  mutually  orthogonal  and  all  of  unit  length) 
x0,  y 0,  z 0.  These  are  rotated  by  ip  about  z0  to  give  ccl5  y1,z1.  These  new  are  now  rotated 
by  9  about  yl  to  give  x2,y2,z2.  Finally,  these  latest  are  rotated  by  cp  about  x2  to  give 
x3,y3,z3: 


IQi)  X"Vl'Zl  SOT  %’y2'*2  V?  X3'V3'Z 3'  (<U9) 


The  question  is:  given  x0,y0,z0  and  x3,y3,z3,  what  are  ip,9,(p7  These  angles  are  as 
follows.  The  angle  ip  is  given  by  projecting  x3  onto  the  x0y0  plane,  and  is  the  angle  of 
this  projection  from  x0  in  the  y0  direction: 


sin  ip 


x3-Vo 

^(x3-x0)2  +  (x3-y0)2  ’ 


cos  ip 


_ XYX0 _ 

V(X3-Xo)2  +  (X3-Vo)2 


(4.20) 


Next,  9  is  the  angle  between  x3  and  the  projection  just  mentioned: 


sin0  =  —x3-Zq  ,  cos9  =  ^(x3-x0)2 +  (x3-yQ)2  .  (4.21) 

Finally,  cp  is  the  angle  from  y2  to  y3  in  the  z2  direction: 

sin  cp  =  y3-z2,  cos  <p  =  y3-y2.  (4.22) 

Matlab  code  to  calculate  the  first  two  of  these  angles  would  be 


psi  =  atan2(  dot(x3,y0),  dot(x3,x0)  ); 

theta  =  atan2(  -dot(x3,z0),  sqrt(  (dot(x3,x0))~2  +  (dot(x3,y0))~2)  ); 

For  the  last  angle  cp,  we  need  y2  and  z2.  These  are  given  by 

y2  =  y1  =  Rzo  (ip)y0  ,  z2  =  Ryi  ( 9)z1  =  RV2  (9)z0  .  (4.23) 

Appropriate  Matlab  code  is  then 
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phi  =  atan2(  dot(y3,z2),  dot(y3,y2)  ); 

The  sequence  of  steps  to  convert  to  the  six  DIS  numbers  from  a  lat-long-height  and 
heading-pitch-roll  is  as  follows: 

1.  Use  (2.2)  to  convert  lat-long-height  (a,ui,h)  to  ( X,Y,Z ). 

2.  Find  the  local  NED  axes  as  we  did  in  Section  4.1.  That  is,  start  with  X,  Y,  Z  in  the 
ECEF  frame  and  rotate,  first  through  the  longitude  and  then  through  the  latitude, 
to  give  —D,E,N,  respectively. 

3.  Rotate  the  AT",  E,  D  vectors,  first  by  the  heading  around  D,  then  by  the  pitch  around 
the  rotated  E,  and  finally  by  the  roll  around  the  latest  rotated  AT".  The  resulting 
vectors  are,  respectively,  the  aircraft’s  local  x,y,z  axes  (all  in  the  ECEF  frame). 

4.  The  angles  are  calculated  from  (4.20,  4.21,  4.22),  where  we  make  the  identi¬ 

fications 


x0  —  X  ,  y0  —  Y  ,  z0  —  Z  , 

x3  =  x,  y3  =  y,  Z3  =  z.  (4.24) 

Example  1:  An  aircraft  is  flying  10,000m  above  Adelaide,  heading  south-east  (with  no 
wind  component:  the  heading  is  taken  as  the  direction  in  which  the  nose  points),  climbing 
at  a  20°  pitch,  and  holding  a  30°  roll.  What  are  the  two  sets  of  triplets  that  the  DIS 
standard  requires  to  specify  the  aircraft’s  location  and  orientation? 

The  location  is  found  easily,  by  converting  the  aircraft’s  lat-long-height  at  Adelaide 
to  XYZ  using  (2.2): 


(X,Y,Z)  =  (-3.93,  3.48,  -3.63)xl06m.  (4.25) 

The  orientation  vectors  x,  y,  z  are  found  by  first  calculating  the  local  NED  axes  as 
vectors  in  the  ECEF,  and  then  rotating  these  by  the  relevant  heading,  pitch  and  roll. 
Some  of  this  calculation  has  been  done  already  in  Section  4.1;  the  local  NED  axis  vectors 
are  given  there  in  (4.8).  Using  these,  place  the  aircraft  into  the  NED  frame  so  that  initially 
its  x,  y,  z  axes  (i.e.  x,  y,  z  vectors)  coincide  with  the  AT",  E,  D  axes  respectively,  and  now 
rotate  it  as  needed.  So  first  set 

x0  =  N,  y0  =  E,  z0  =  D.  (4.26) 

Now  rotate  x0,y0,z0  each  by  135°  around  z0  (implementing  the  south-east  heading), 
calling  the  result  x1,y1,z1.  Then  rotate  each  of  these  by  20°  around  y1  (implementing  the 
pitch)  to  produce  x2,  y2,  z2.  Finally  rotate  each  of  these  by  30°  around  x2  (implementing 
the  roll)  to  produce  x,y,z.  Each  of  these  rotations  is  easy  to  work  out  using  the  basic 
formula  (3.7).  The  result  is 


'-0.366' 

0.928' 

0.065' 

X  = 

-0.564 

,  y  = 

-0.165 

,  *  = 

-0.809 

-0.741 

-0.333 

0.584 
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Now  calculate  tj},9,(j>  from  (4.20,  4.21,  4.22).  To  do  this,  set 


II 

o 

i 

O 

_ i 

II 

o 

e* 

i 

i- 1  O 

_ i 

II 

o 

fc* 

1 - 

O  O 

i _ 

0 

0 

l 

x3  =  x,  y3  =  y ,  z3  =  2:-  (4.28) 

The  resulting  angles  are 

^  =  -123.0°,  9  =  47.8°  ,  (/)  =  - 29.7°.  (4.29) 

These,  together  with  (X,  Y,  Z)  given  in  (4.25),  are  the  six  numbers  that  the  DIS  standard 
uses  to  encode  the  aircraft’s  position  and  orientation. 

Example  2:  Convert  back  from  X,Y,  Z,ip,9,(j>  to  lat-long-height/heading-pitch-roll. 

Again  the  location  is  found  easily  by  applying  (2.3)-(2.6).  These  equations  return  the 
latitude  a  =  —34.9°,  longitude  oj  =  138.5°,  and  height  10,000m,  as  expected. 

To  calculate  the  heading,  pitch,  and  roll,  realise  that  these  are  just  the  'i/j,  9,  cj),  respec¬ 
tively,  that  are  returned  by  step  4  above  (4.24)  when  we  make  the  following  identifications: 

x0  =  N ,  y0  =  E  ,  z0  =  D  , 

x3  =  x,  y3  =  y,  Z3  =  Z.  (4.30) 

So  we  need  to  calculate  N ,  E,  D,  which  can  be  done  as  in  Section  4.1,  since  we  have  just 

calculated  the  latitude  and  longitude  of  the  aircraft.  These  are  given  in  (4.8). 

Next,  the  x,y,z  are  found  by  applying  the  three  Euler  rotations  of  tf) ,9,(f>  in  the 
appropriate  order,  i.e.  (4.19),  to  the  ECEF’s  basis  vectors  X,Y,Z.  In  other  words,  set 
x0,  y 0,  z 0  in  (4.19)  to  X,  Y,  Z.  The  three  resulting  vectors  x3,  y3,  z3  will  be  the  required 
x,  y,  z. 

Finally,  make  the  identifications  of  (4.30)  to  calculate  il’,9,(f)  using  (4.20,  4.21,  4.22), 
where  the  heading  is  i/j,  the  pitch  is  9 ,  and  the  roll  is  4>.  As  expected,  we  obtain 

heading  =  135°  ,  pitch  =  20°  ,  roll  =  30°  .  (4-31) 


5  Concluding  Remarks 

All  of  the  commonly-needed  techniques  in  three-dimensional  rotations  rely  on  our  ability 
to  do  two  basic  procedures: 

1.  Rotate  a  vector  about  any  axis  whatsoever,  and 

2.  calculate  components  using  dot  products. 
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These  calculations  need  to  be  done  in  just  one  frame.  The  usual  one  is  the  ECEF,  which 
is  good  for  all  of  the  usual  motions  we  associate  with  moving  around  Earth.  (If  we  need 
to  deal  with  satellites,  a  nonrotating  frame  is  needed,  such  as  one  fixed  to  the  stars.) 
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Rotating  about  arbitrary  axes  and  taking  dot  products  within  the  ECEF  is  all  that  has 
been  done  for  the  examples  in  this  report. 

Rotations  themselves  can  be  carried  out  using  matrices  or  quaternions.  Both  of  these 
are  simply  ways  of  writing  down  the  three  numbers  needed  to  encode  a  rotation  (an 
angle,  and  two  numbers  for  the  axis).  The  extra  redundancy  that  matrices  have  over 
quaternions  allows  for  simpler  algebra  when  using  matrices  to  analyse  rotation  on  paper, 
but  quaternions,  being  essentially  pared-down  versions  of  the  rotation  matrix,  can  be  easier 
to  prevent  from  degrading  numerically  inside  a  computer  routine  (but  see  the  comment  at 
the  end  of  Sect.  3.4).  Performance  questions  like  this  are  critical  in  that  three-dimensional 
models  can  have  10,000+  points  to  rotate  at  e.g.  30  frames  per  second. 

When  solving  a  problem  involving  three-dimensional  rotations,  a  good  approach  is  to 
break  the  scenario  down  into  its  building  block  rotations  and  to  handle  each  one  separately, 
all  within  the  ECEF.  This  step-by-step  approach  lends  itself  well  to  being  modularised  in 
software,  and  is  what  we  have  used  repeatedly  in  this  report. 
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Appendix  A  Combining  Two  Rotations 

Earlier  we  said  that  two  rotations  will  always  combine  to  give  a  new  rotation.  This  is  of 
course  equivalent  to  multiplying  the  two  rotation  matrices  or  quaternions.  But  what  is 
the  axis  and  angle  of  the  new  rotation?  For  example,  suppose  we  rotate  a  vector  twice: 

1.  First  rotate  it  by  90°  around  the  x-axis  (call  the  matrix  Rx  and  the  quaternion  Qx). 

2.  Then  rotate  it  by  90°  around  the  y-axis  (call  the  matrix  Ry  and  the  quaternion  Qy ). 

What  is  the  corresponding  axis  and  angle  of  the  combined  rotation?  We  do  this  first  using 
matrices  and  then  with  quaternions.  It  will  be  quite  evident  from  the  answer  that  two 
rotations  generally  don’t  combine  to  give  an  obvious  resulting  angle  and  axis. 


Using  Matrices 


The  combined  matrix  rotation  is  Ry  Rx : 


T 

T 

0 

O' 

1.  Rx  :  8 

=  90°,  n  = 

0 

,  so  Rx  = 

0 

0 

-1 

0 

_0 

1 

0_ 

'O' 

0 

0 

1 

2.  Ry  :  8 

=  90°,  n  = 

1 

,  so  Ry  = 

0 

1 

0 

0 

-1 

0 

0 

3.  Ry  Rx  = 


0  1  0 
0  0-1 
-10  0 


To  find  the  single  angle  and  axis  of  rotation  that  Ry  Rx  is  equivalent  to,  we  could  compare 
each  element  of  Ry  Rx  with  the  general  expression  (3.7),  but  this  is  arduous.  Much  easier 
is  to  use  (3.7)  to  write  down  two  properties  of  the  general  rotation  matrix,  using  its  trace 
and  its  transpose: 


tr  R  =  1  +  2  cos  8 

R-Rl  =  2  sin  0nx. 


(Al) 


In  that  case,  the  combined  rotation  is  around  some  unit  vector  n  through  angle  8,  where 


l  +  2cos0  =  O  ,  2  sin  8nx 


0  1  1 

-1  0  -1 

-1  1  0 


(A2) 


There  appear  to  be  two  solutions  for  n  and  8  here,  but  they  both  represent  the  same 
rotation,  so  that  we  can  always  choose  the  unique  positive  value  of  8  together  with  its 
corresponding  n.  In  this  case: 

8  =  120°  ,  n  =  (  1, 1,  — 1)/V3 .  (A3) 
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Using  Quaternions 

The  necessary  quaternions  are 

Qx  =  (cos 45°,  sin 45° (1,0,0))  =  -)=(1,  1,0,0) , 
Qy  =  (cos  45°,  sin  45°  (0,1,0))  =  -^=(1,0, 1,0) . 


The  combined  quaternion  rotation  is  Qv  Qx : 


(A4) 


Qy  Qx  =  |(1,0, 1,0)  (i,  1,0,0)  =  |(i,  l,  l, -l) 

(l  (  cno  (i,i,-i)  .  enQ1 

=  7=-j  =  ^cos60,— ^smOO  ), 

which  is  a  rotation  of  120°  around  (1, 1,  —  1)/ a/3,  just  as  we  found  in  (A3) 


(A5) 
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Appendix  B  xyz  Eulers  to  Angle— Axis  and  Vice 

Versa 

In  this  appendix  we  explain  how  to  switch  between  an  xyz  Euler  angle  representation  of  a 
rotation  and  the  angle-axis  representation  characterised  by  a  single  rotation  matrix  Rn(9) 
or  quaternion  Qn(9).  This  fixed-axis  Euler  representation  differs  to  the  DIS  standard 
discussed  in  Section  4.3  because  in  this  appendix,  the  axes  rotated  around  are  fixed  in 
space.  Even  so,  the  fixed-axis  convention  is  also  very  commonly  used. 


xyz  Eulers  to  Angle— Axis 


Suppose  a  rotation  is  built  of  three  Euler  rotations:  first  by  an  angle  a  around  the  x-axis, 
then  by  (5  around  y,  then  by  7  around  z.  These  join  to  give  a  single  rotation  of  9  around  n, 
given  by 

Rn{e)  =  Rz(1)Ry(P)Rx{a ),  (Bl) 

where  the  matrix  Rn{9)  encodes  the  whole  rotation.  If  it’s  necessary  to  know  n  and  9 
explicitly,  they  can  be  found  by  applying  (Al). 

If  the  rotations  are  performed  around  new  axes,  then  the  calculation  of  the  rotation 
matrix  is  only  a  little  more  difficult.  For  example,  suppose  the  second  rotation  is  around 
the  new  y-axis  (called  y'),  while  the  third  is  around  the  newest  x,  called  x': 


Rn{9)  =  RX’{ 7)  Ry’(/3)  Rx{a) . 


(B2) 


The  first  rotation  is  done  as  usual:  Rx{a)  rotates  by  a  around  (1, 0, 0).  The  second  rotation 
is  about  the  new  y-axis,  whose  vector  is 


nx  =  Rx(a ) 


0 

1 

0 


(B3) 


The  last  rotation,  by  7,  needs  to  be  done  around  the  new  x-axis.  This  axis  is  represented 
by  the  vector 


n2  =  Rni(P) 


1 

0 

0 


(B4) 


so  that  the  final  rotation  is 


Rn{9)  =  Rn2 (7)  Rni  (P)  Rx{a) . 


(B5) 


These  calculations  can  of  course  also  be  done  using  quaternions,  as  detailed  in  Appendix  A. 
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xyz  Angle— Axis  to  Eulers 


If  the  xyz  rotation  order  of  (Bl)  is  used,  a  single  rotation  matrix  R  can  be  converted  to 
three  rotations  using  the  Euler  angles  a,/3, 7,  where  s  can  be  chosen  as  either  +1  or  —1: 


sin  ct  = 

SR32 

cos  a  = 

SR33 

yi  -  Am  ' 

Vl-^31 

sin  (3  = 

—  R'31  , 

cos  (3  = 

s\/1  ~  ^31 

sin  7  = 

si?2l 

cosy  = 

sRn 

yi  -  ' 

Vi-^1 

(B6) 


The  value  of  s  is  immaterial  in  the  sense  that  either  choice  will  give  three  angles  that 
have  the  same  effect  as  the  original  rotation.  Care  is  needed  to  avoid  oversimplifying  these 
equations:  we  must  remember  that  (B6)  does  not  imply  that  e.g.  a  =  tan”1  (R32/R33), 
since  the  tan”1  function  only  gives  angles  in  the  range  —  90°  — >  +90°.  Nor  does  it  imply 
that  P  must  equal  —  sin”1  R31,  because  of  a  similar  restriction  on  the  sin”1  function. 


These  equations  can  be  implemented  in  Matlab  using  the  code: 


alpha 

beta 

gamma 


atan2(s*R(3,2) , 
atan2(-R(3, 1) , 
atan2 (s*R(2 , 1) , 


s*R(3,3)) ; 

s*sqrt (1-R(3, 1) ~2) ) ; 
s*R(l,l)) ; 
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Appendix  C  Confusing  Euler  Angle  Orientation 

with  Incremental  Rotation 

As  we  saw  in  Sect.  3.2,  an  object’s  orientation  can  be  specified  by  giving  three  Euler  angles: 
angles  through  which  the  object  can  be  turned  in  some  preset  order  from  a  base  position, 
which  will  result  in  the  required  orientation.  This  is  useful  and  unambiguous,  but  it  can 
lead  to  some  confusion  when  implemented  in  an  overly  rigid  way. 

The  problem  can  be  shown  as  follows.  Suppose  we  have  written  computer  code  that 
draws  some  object,  and  supplies  three  graphical  sliders  that  allow  us  to  specify  Euler  angles 
to  change  its  orientation.  We  can  alter  any  of  the  angles  at  any  time,  and  the  software  will 
read  the  current  values  of  the  three  angles,  interpret  them  as  Euler  angles,  and  use  them 
to  rotate  a  preset  base  orientation  of  the  object,  first  around  the  x-axis,  then  around  y, 
and  finally  around  z.  The  result  is  then  displayed  on  the  computer  screen. 

When  we  start  the  programme,  this  graphical  interface  is  always  drawn  with  the  three 
sliders  each  set  to  zero,  and  the  object  in  its  base  position.  If  we  alter  the  x-slider  to 
read  10°  and  leave  the  other  two  sliders  at  zero,  then  the  three  angles  are  read  by  the  soft¬ 
ware,  and  a  rotation  of  Rz( 0)  Ry{ 0)  i?x(10°)  is  applied.  The  object  duly  rotates  around  x 
by  10°.  If  we  increase  the  x-slider  to  15°,  the  software  reads  all  the  sliders  again  and 
applies  a  rotation  of  Rz(0)  Ry(0)  Rx(15°)  to  the  base  orientation.  The  effect  is  that  the 
object  is  seen  to  rotate  around  x  through  a  further  5°.  We  then  set  the  x-slider  back  to 
zero,  the  three  angles  are  again  read,  and  the  result  is  that  the  base  orientation  is  rotated 
by  Rz(0)  Ry(0)  Rx(0),  or  nothing  at  all.  So  changing  the  x-slider  by  itself  rotates  the 
object  about  the  x-axis,  which  is  perfectly  reasonable  and  expected. 

Similarly  if  we  change  only  the  y-slider,  or  only  the  z-slider,  the  same  sort  of  intuitive 
effects  happen.  Setting  the  y-slider  to  10°  with  the  others  at  zero  means  that  the  software 
reads  the  three  angles  and  applies  a  rotation  of  Rz( 0)  ^(10°)  -Rx(0)  to  the  base  orientation, 
so  that  the  object  rotates  around  y  by  10°,  and  similarly  it  will  revert  to  the  base  position 
if  we  set  the  y-slider  back  to  zero. 

Suppose  we  now  change  two  of  the  angles;  we’ll  ignore  the  z-slider  as  it’s  not  needed  to 
make  the  following  point.  Start  with  all  sliders  at  zero  so  that  the  object  is  drawn  in  the 
base  orientation.  Then  set  the  x-slider  to  10°.  The  rotation  is  Ry{ 0)  i2a,(10°),  so  the  object 
rotates  around  x  by  10°.  Now  set  the  y-slider  to  10°.  The  new  rotation  is  -R^IO0)  i2x(10°) 
from  the  base  orientation,  with  the  result  that  the  object  is  first  rotated  around  x  by  10° 
and  then  around  y  by  10°.  So  far,  the  behaviour  matches  our  intuition. 

Now  increase  the  x-slider  to  read  15°.  The  rotation  is  now  ^(10°)  Rx{  15°)  from  the 
base  orientation,  so  that  the  object  is  first  rotated  around  x  by  15°  and  then  around  y 
by  10°.  The  procedure  is  shown  as  follows: 


x-slider 

setting 

y-slider 

setting 

Resulting  rotation 
of  base  orientation 

0 

0 

Ry{ 0)  ik(0) 

10° 

0 

Ry{ 0)  Rx(  10°) 

10° 

10° 

Ry( 10°)  Rx( 10°) 

15° 

10° 

Ry( 10°)  Rx( 15°) 
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But  the  effect  of  this  latest  rotation  is  not  what  we  might  have  expected.  It  is  reasonable 
to  believe,  based  on  trying  each  slider  in  turn  on  its  own,  that  each  slider  rotates  about  its 
respective  axis.  When  we  increased  the  x-slider  from  10°  to  15°,  most  users  would  expect 
the  object  to  rotate  by  a  further  5°  about  x.  But  this  is  not  what  happens.  The  software 
does  what  it  was  designed  to  do:  it  rotates  the  base  orientation  using 

i2y(10°)  ik(15°).  (Cl) 

The  rotation  suggested  by  our  intuition  is  an  increment  of  5°  about  x  on  the  last  orientation — 
which  was  represented  by  72^(10°)  i?x(10°)  acting  on  the  base  orientation — or 

Rx(§°)  -Rj/(10°)  i?x(10°) .  (C2) 

Expressions  (Cl)  and  (C2)  are  not  the  same.  In  fact  (Cl)  could  have  been  written 

i*v(10°)  Rx( 5°)  74(10°) ,  (C3) 

which  differs  from  (C2)  in  that  two  rotations  are  swapped.  Because  rotation  is  noncom- 
mutative,  these  two  different  matrices  (C2)  and  (C3)  produce  different  orientations  when 
multiplying  the  base  orientation.  To  reiterate,  our  intuition  expected  that  by  making  three 
slider  adjustments,  first  incrementing  the  x-slider  (from  zero)  by  10°,  then  incrementing 
the  y-slider  (from  zero)  by  10°,  then  incrementing  the  x-slider  by  5°,  that  (C2)  would 
result.  In  fact,  (Cl)  or  equivalently  (C3)  resulted.  That  might  be  unintuitive,  but  the 
mathematics  and  the  computer  software  have  done  exactly  what  was  asked  of  them. 


Our  confusion  over  the  effects  of  the  sliders  seems  to  reach  its  worst  if  the  y-slider  is 
first  set  to  rotate  through  —90°,  and  then  the  x-  and  2-sliders  are  changed  arbitrarily. 
Because  the  following  identity  holds  for  all  angles  ct,  (3 , 


Rz(a)  Ry(- 90°)  Rx{(3) 


0  —  sin(a  +  /3)  —  cos(a  +  (3 ) 

0  cos(a  +  (3)  —  sin(a  +  (3)  , 

1  0  0 


(C4) 


the  total  rotation — being  a  function  only  of  the  sum  of  the  x-  and  2-slider  angles — does 
not  distinguish  between  the  x-  and  2-sliders.  Thus  the  x-  and  2-sliders  have  the  same 
effect,  making  it  appear  that  we  have  lost  one  rotation  axis  somewhere,  as  if,  to  use  an 
oft-quoted  phrase,  “the  x-axis  has  been  rotated  onto  the  2- axis”  (!)  But  axes  are  not  set 
on  hinges;  there  is  no  such  thing  as  a  flexible  axis,  and  axes  cannot  be  rotated  onto  each 
other.  The  software  is  doing  exactly  what  it  was  designed  to  do.  Whether  our  intuition 
would  prefer  a  different  behaviour  is  another  matter  entirely. 

This  apparent  loss  of  one  axis  is  erroneously  labelled  gimbal  lock  by  some  graphics 
programmers,  who  confuse  it  with  the  real  numerical  Euler  angle  problem  described  earlier 
on  page  7.  But  there  is  no  such  problem  with  Euler  angles  in  software  rotations  such  as 
described  above,  and  this  use  of  the  term  is  a  misnomer.  Unfortunately,  historically  these 
misinterpretations  went  hand  in  hand  with  using  matrices  to  perform  rotations,  with  the 
result  that  matrices  and  Euler  angles  are  sometimes  seen  as  not  up  to  the  task.  This  is 
definitely  wrong;  both  matrices  and  Euler  angles  can  handle  any  rotation  asked  of  them. 
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Restoring  an  Intuitive  Feel  to  the  Sliders  If  we  wish  the  sliders  to  act  in  an  intuitive 
way,  then  each  time  we  adjust  any  slider,  we  must  not  have  the  software  read  the  three 
slider  values  and  rotate  the  base  orientation  through  them  all  in  some  preset  order.  Rather, 
if  we  adjust  the  x-slider  by  2°,  the  software  should  read  the  values,  note  what  has  changed, 
and  rotate  the  current  orientation  about  the  x-axis  through  the  increment  of  2°.  Not 
only  does  this  make  the  sliders  and  rotations  behave  as  our  intuition  suggests,  but  it  also 
makes  for  faster  processing,  since  only  one  rotation  is  ever  performed  per  slider  interaction, 
instead  of  three. 
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Appendix  D  Quaternions  Used  in  Computing: 

SLERP 

Rotation  theory  is  important  when  creating  a  smooth  fly  through  of  a  scene  in  which 
certain  way  points  have  been  specified  along  with  camera  orientations  at  those  points.  A 
camera’s  orientation  can  be  specified  by  a  rotation  matrix  or  quaternion  that  generates 
that  orientation  by  rotating  the  initial  orientation.  So  the  question  arises  as  to  how  best 
to  interpolate  the  matrices  or  quaternions,  in  order  to  generate  rotations  that  act  on  the 
initial  orientation  to  give  visually  acceptable  intermediate  orientations. 

For  example,  if  two  way-point  orientations  specified  by  rotation  matrices  Rq  and  R\ 
are  specified,  then  we  might  suppose  that  linear  combinations  of  these  with  appropriate 
weightings  will  suffice  to  rotate  the  initial  orientation  to  create  intermediate  orientations. 
But  such  is  not  the  case;  the  intermediate  rotation  matrices  are  not  simply  given  by  linear 
interpolations  between  Rq  and  Ri-  For  example  in  the  simpler  case  of  rotations  in  the 
xy  plane,  if  the  rotation  matrix  for  30°  is  -R2(30°)  using  (3.4),  and  that  for  40°  is  i?2(40°), 
then  the  matrix 

0.9R2(30°)  +  0.1R2(40o)  (Dl) 

is  not  equal  to  i?2(31°);  being  not  orthogonal,  it’s  not  a  rotation  matrix  at  all.  Even 
when  orthogonalised  using  (3.13),  what  results  still  does  not  equal  i?2(31°).  So  a  more 
sophisticated  matrix  interpolation  is  needed. 

Whatever  results  are  applicable  to  matrices,  quaternions  seem  to  be  easier  to  use; 
and  certainly  most  of  the  research  into  interpolating  orientations  seems  to  revolve  around 
quaternions.  This  seems  to  be  mostly  due  to  the  fact  that  because  of  its  numerically 
intensive  nature,  rotation  research  tends  to  be  mainly  the  domain  of  computer  animation 
programmers;  but  unfortunately  once  various  quaternion  routines  had  been  worked  out 
and  coded,  the  code  got  copied,  understood,  and  used  extensively;  and  quaternions  became 
the  de  facto  way  to  solve  problems  of  this  sort.  Nonetheless,  as  we’ll  show  in  this  appendix, 
quaternions  are  very  easily  interpolated  to  give  more  acceptable  interpolation  results  than 
are  traditionally  obtained  using  matrices. 

The  fact  that  quaternions  have  a  unit  length  given  by  a  simple  sum  of  squares  of  their 
components  means  they  can  be  treated  as  vectors  in  Euclidean  four  dimensions,  so  are 
able  to  be  visualised  as  points  on  the  surface  of  a  “sphere”,  as  shown  in  Figure  Dl.  In  this 
figure,  the  initial  orientation  is  represented  by  the  quaternion  (1,0,  0,0),  because  this  is 
the  identity  quaternion  that  doesn’t  rotate  at  all,  shown  by  applying  (3.17): 

(0,  x')  =  (1,0, 0, 0)  (0,  *)  (1,  0,  0,  0)  =  (0,  x) .  (D2) 

We  wish  to  generate  a  whole  series  of  quaternions,  each  of  which  multiplies  the  initial 
orientation  to  produce  some  intermediate  one.  We  must  pass  through  the  way-point 
quaternion  (0.9,  0.1,  0.1, 0.4)  and  finish  at  the  final  quaternion  (0.7, 0.6,  0.2, 0.3).  The  very 
simplest  way  is  to  make  a  linear  interpolation,  in  angle,  between  way  points  on  the  sphere’s 
surface,  and  so  is  called  spherical  linear  interpolation ,  or  slerp.  To  see  how  slerp  is  used 
to  interpolate  between  two  quaternions,  draw  them  as  vectors  q1  and  q2  lying  in  a  plane 
in  Figure  D2.  (Previously  we  used  e.g.  Q  to  denote  a  quaternion.  Here  we  use  bold  face 
vector  notation  to  reinforce  our  treatment  of  quaternions  as  vectors,  since  that’s  what 
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Figure  Dl:  Schematic  of  a  four- dimensional  sphere,  each  surface  point  of 
which  represents  a  quaternion.  Any  orientation  of  some  system  is  given  by 
a  quaternion  that  rotates  the  initial  orientation.  Given  way-point  orienta¬ 
tions,  the  task  is  to  find  a  path  connecting  them  that  then  picks  out  the  set 
of  interpolating  quaternions  that  each  generate  an  acceptable  intermediate 
orientation  of  the  system 


picturing  them  as  points  on  the  surface  of  a  four- dimensional  sphere  allows  us  to  do.)  The 
angle  9  between  them  is  given  by  the  usual  dot  product: 

cos  9  =  Ql  92  =  q1q2  ,  (D3) 

|9il  1921 

where  the  dot  product  is  the  usual  sum  of  pairwise  products  of  components.  Use  a  param¬ 
eter  t  running  from  zero  to  one,  values  of  which  generate  intermediate  quaternions  such 
as  the  q  shown  in  Figure  D2.  Straightforward  two-dimensional  vector  analysis  gives 


sin(l  —  t)9  sin  t9 

9  =  ZTTTn  9l  +  :  7T  92 


sin# 


sin# 


(D4) 


This  is  the  standard  SLERP  formula  for  generating  intermediate  quaternions.  In  Figure  Dl 
we  first  interpolate  for  various  values  of  t  between  the  initial  orientation  and  the  way  point, 
and  then  do  the  same  for  the  way  point  and  the  final  orientation.  Remember  that  each 
quaternion  generated  does  not  multiply  the  orientation  just  produced;  rather,  it  multiplies 
the  initial  orientation. 


Figure  D2:  The  quaternion  q  is  interpolated  between  quaternions  g1 


and  q2 
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SLERP  is  popular  among  graphics  programmers,  but  there  is  some  recognition  that 
this  is  partly  by  default  because  of  the  way  routines  get  shared  around.  SLERP  is  not  the 
only  quaternion  tool  used,  and  it  has  even  been  argued  as  generally  not  desirable  to  use, 
partly  because  of  its  high  processing  demands,  and  partly  because  other  algorithms  (such 
as  normalised  quaternion  linear  interpolation)  bring  other  benefits  to  the  programmer  in 
terms  of  the  “look  and  feel”  of  the  resulting  rotations.  Nevertheless,  the  above  description  of 
SLERP  shows  the  general  idea  of  the  current  approach  to  interpolation  within  the  computer 
graphics  community. 
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